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A single-level quantum dot with Coulomb repulsion attached to two superconducting leads is 
studied via the perturbation expansion in the interaction strength. We use the Nambu formalism and 
the standard many-body diagrammatic representation of the impurity Green functions to formulate 
the Matsubara self-consistent perturbation expansion. We show that at zero temperature second 
order of the expansion in its spin-symmetric version yields a nearly perfect agreement with the 
numerically exact calculations for the position of the 0 — tt phase boundary at which the Andreev 
bound states reach the Fermi energy as well as for the values of single-particle quantities in the 
0-phase. We present results for phase diagrams, level occupation, induced local superconducting 
gap, Josephson current, and energy of the Andreev bound states with the precision surpassing any 
(semi)analytical approaches employed thus far. 


INTRODUCTION 

Nanostructures attached to leads with specific prop¬ 
erties display interesting and important quantum effects 
at low temperatures. Much attention, both from exper¬ 
imentalists m and theorists 0, has been paid in recent 
years to a quantum dot with well separated energy lev¬ 
els attached to BCS superconductors. In particular, the 
behavior of the supercurrent (Josephson current) that 
can flow through the impurity in equilibrium without 
any external voltage bias between two superconducting 
leads was in the center of interest nn. The Joseph¬ 
son current through quantum dots with tangible on-dot 
Coulomb repulsion can induce a transition signalled by 
the sign reversal of the supercurrent observed experimen¬ 
tally |6l-[l^. 

This so called 0 — tt transition is induced by the un¬ 
derlying impurity quantum phase transition (QPT) re¬ 
lated to the crossing of lowest many-body eigenstates of 
the system from a spin-singlet ground state with positive 
supercurrent (0-phase) to a spin-doublet state with neg¬ 
ative supercurrent (7r-phase) |ll| . In single-particle spec¬ 
tral properties this transition is associated with crossing 
of the Andreev bound states (ABS) at the Fermi energy 
as has also been experimentally observed dllil. Con¬ 
tinuous vanishing of the ABS energies at the transition 
is a direct consequence of crossing of many-body eigen¬ 
states and may serve as an important consistency 
check of proposed theories. The latter cover by now a 
broad scope of techniques ranging from numerically exact 
(and computationally expensive) numerical renormaliza¬ 
tion group (NRG) (l4 Il5j suitable for zero-temperature 
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and finite-temperature quantum Monte Carlo (TH, [T^ to 
(semi)analytical methods based on expansion around the 
atomic limit [Tsl - I^ . mean-field theory (2il - l^ . or for¬ 
malisms specialized on the strongly correlated systems 
such as slave-particl es and functional renormal¬ 

ization group (fRG) j^. 

However, despite of the versatility of these approaches, 
there still remain vast regions of the parameter space 
with direct experimental relevance (A < T < {7, see, 
e.g.. Ref. @) where most of the above approaches can¬ 
not be applied and one has to resort either to overly 
heavy numerical methods (NRG or QMC) or to conceptu¬ 
ally flawed spin-symmetry-broken mean-field approach. 
The latter approach is not excessively elaborate and of¬ 
ten gives quantitatively acceptable results jl^, although 
at the expense of breaking the spin symmetry of exact 
solution. In particular, spin-polarized mean-field solu¬ 
tions even after the symmetrization described in Ref. 
still exhibit at the transition unphysical discontinuities in 
the ABS energies jU and in finite-temperature supercur- 
rents 

The aim of this paper is to provide a conceptually clean 
and computationally inexpensive generic formalism for 
addressing the 0 — tt transition in that widespread regime 
without strong-correlations (i.e., without the Kondo ef¬ 
fect). We show that a resummed perturbation theory 
(PT) incorporating the second-order dynamical correc¬ 
tions to the spin-symmetric mean-field (Hartree-Fock) so¬ 
lution yields at zero temperature a nearly perfect descrip¬ 
tion of the 0-phase including the position of the phase 
boundary in a wide parameter range outside of strong 
correlations. The precision of this solution is unprece¬ 
dented by any so far employed (semi) analytical methods 
including fRG. On the other hand, the solution developed 
from the non-interacting limit breaks down at the phase 
boundary and any perturbative description of the tt- 
phase and, consequently, also finite temperatures which 
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mix 0 and tt solutions, remains elusive. Although the 
second-order PT has been applied to this problem previ¬ 
ously in Ref. [U Sec. V] and, especially, in a. (otherwise 
unpublished) part of Meng’s master thesis |^L Ch. 4], 
these studies were limited to the particle-hole symmet¬ 
ric case only (in Ref. [Isj in just two limits A/P <C 1 
and A/r ^ 1) and did not use crossing of the ABS as 
the definition of the boundary of the 0-phase. Instead, 
they defined the 0 — tt transition by equalling the ap¬ 
proximated Kondo temperature with the superconduct¬ 
ing gap, namely A = r/(l — i9I](0)/5a;), which however 
holds only qualitatively. The generic character of the PT 
method and the proper definition of the 0 — tt bound¬ 
ary in the Green-function formalism have thus remained 
unnoticed. 


RESULTS 

A single impurity Anderson model is used to simulate 
the quantum dot with well-separated energy levels con¬ 
nected to the superconducting leads in the experimental 
setup [mil [Tr [ 2 ^ . The Hamiltonian of the system 
consisting of a single impurity with the level energy e 
and local Coulomb repulsion U attached to two super¬ 
conductors reads 

R = e ^ dU. + Ud\d/^d^+ ^ (1) 

cr=t 54 ' S—R,L 

where the BCS Hamiltonian of the leads is 

"^fead = +H.C.) , 

kcT k 

( 2 ) 

with s = L,R denoting the left/right lead, respectively. 
Finally, the hybridization term between the impurity and 
the contacts is given by 

H?^ = -i.E(4X+H.c.) . (3) 

kcr 

The individual degrees of freedom of the leads are 
unimportant for the studied problem and are generally 
integrated out, leaving us with only the active variables 
and functions on the impurity. Due to the proximity 
effect there are locally induced superconducting correla¬ 
tions on the impurity and the most direct way to handle 
them is via the Nambu spinor representation of the local 
fermionic operators in which the one-electron impurity 
(imaginary time/Matsubara) Green function (GF) is a 
2x2 matrix 


Gair - t') = 


Gair-r'), 

Ga{T-T'), G-a{T-T') 

/ (TK(r)4(r')]) , (TK(r)d_4r')]) ' 

V(T[d4(T)4(r')]) , {T[dl^{T)d_^{T')])^ 

(4) 


where the bar denotes the hole function. 

The impurity GF can be exactly found for an impurity 
without onsite interaction (U = 0) by method analogous 
to Appendix A of Ref. [l^ . When assuming identical left 
and right superconducting gaps A^ = A/j = A as well as 
tunnel couplings = tn = t it can be written in terms 
of Matsubara frequencies a;„ = {2n + l)TT/f3 as (e = h = 1 
throughout the paper; we also skip the spin index as we 
only consider spin-symmetric solutions) 


- + s{iuJn)] - e , A$s(iw„) 

° V A$s(iw„) , iuin[l + s{iuJn)] + £ 

(5) 

where sUuJn) = —r=- is the hybridization self-energy 

due to the coupling of the impurity to the superconduct¬ 
ing leads. We have denoted by F = 27rt^po th® normal- 
state tunnel coupling magnitude (po being the normal- 
state density of states of lead electrons at the Fermi en¬ 
ergy) and A<i, = Acos(<I>/2) with 4) = 4)^ — 4>i{ being 
the difference between the phases of the left and right 
superconducting leads. 

The impact of the Goulomb repulsion U > 0 on the 
Green function is included in the interaction self-energy 

matrix Tt{iu!n) = ^ propa¬ 

gator in the spin-symmetric situation is determined by 
the Dyson equation G~^{iu}n) = GQ^{iu}n) — T,{iu!n)- 
The symmetry relations for the Green function equa¬ 
tion (III) reformulated in the Matsubara representation as 
GaiioJn) = -Ga{-iuJn) and Ga{i^_n) = Ga{-i^n) imply 
the same for the self-energies, i.e. Ecr(ia;„) = —Eo.(—iw^) 
and Sa{iuJn) = Sa{—iuJn)- Therefore, the Green function 
explicitly reads 



G{iuJn) = - 


1 


D{iuJn) 

ioJnli + s{iuJn)] -I- e -I- S(-ia;„), 
-A$s(iw„) -|-5(-ja;„). 


iUn 


-A$s(ia;„) -I- S{iu}n) 

[1 -I- s(iw„)] - £ - E(ia;„) 
( 6 ) 


The negative determinant of the inverse 
Green function D(iujn) = — det[G“^(iw„)] = 

[1 + s{iujn)] + [e -|- S(iu;„)] [e -|- S(—ia;„)] -I- 

[A$s(ia;„) — 5(iw„)] [A$s(ia;„) — 5(—ia;„)] deter¬ 

mines via its zeros the existence and positions of the 
ABS. This determinant is real within the gap and can 
go through zero D{ujq) = 0 determining the (real) in-gap 
energies ±wo of the ABS symmetrically placed around 
the Fermi energy (center of the gap). The ABS are 
important for transport of the Gooper pairs through 
the quantum dot and usually provide the dominant 
^ contribution to the dissipation-less Josephson current J 
through the impurity, which can be evaluated at zero 
temperature by an integral of the anomalous Green 
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Figure 1. Diagrammatic representation of the first two orders 
of the perturbation expansion in the Coulomb interaction for 
the normal (top) and anomalous (bottom) parts of the self¬ 
energy. The wavy line represents the Coulomb interaction 
and the lines with single (double) arrow represent the normal 
(anomalous) propagators according to equation (|4)). 


function (see the Methods section) 


4A 



doJn 

27r 


g(iuJn)s(iu;„')e 


—r sin 


$ 


Res(g; -ojo) 
- wg 


duj mjuj) 

7-00 TT - A2 

(7) 


While the first line uses the thermal representation via 
Matsubara frequencies the second one is the analytic con¬ 
tinuation to the real frequencies (spectral representation) 
which allows us to distinguish the direct supercurrent 
through the lower ABS (corresponding to the residue of 
the anomalous impurity Green function at the negative 
ABS frequency) from the tunneling current between the 
continuum band states below the SC gap. 


Spin-symmetric Hartee-Fock approximation 

As the exact expression for this model’s self-energy is 
unknown we resort to the standard Matsubara perturba¬ 
tion theory summing one-particle irreducible diagrams 
for the self-energy. 

The simplest diagrams are the first-order Hartree-Fock 
contributions represented by the first diagrams on the 
r.h.s. of equations in Fig. [TJ Their mathematical equiv¬ 
alents read 

^ G{iuJn) and X! 

n^Z 

( 8 ) 

The HF approximation leads just to a static, 
frequency-independent mean-field self-energy neglecting 
any dynamical correlations caused by particle interac¬ 
tion. Despite of this simplicity and contrary to the com¬ 
mon belief, this approximation yields without any sym¬ 
metry breaking the 0—tt quantum phase transition and we 
thus use it as a convenient and sufficiently simple demon¬ 
stration of the generic features of the perturbation expan¬ 
sion. The Hartree-Fock approximation consists of two 
self-consistent non-linear equations that can be reformu¬ 
lated in terms of auxiliary quantities Ed = e -|- 17 (dj-do-) 
(mean-field energy of the level) and d = F cos($/2) -|- 





e/A 


Figure 2. Phase diagrams in the F — U (a), [/ — <& (b), and 
r — e (c) parameter planes. We compare the phase bound¬ 
aries calculated by numerically exact NRG with various an¬ 
alytical approximations: £RG (only in panel a; data taken 
graphically from Fig. 2 of Ref. |^h spin-symmetric HF, the 
second-order PT/dynamical corrections (DC), and general¬ 
ized atomic limit approximation (GAL) 17^/(1 -I- P/A)^ = 
(2e-b H)^-b dP^ cos2($/2). 


(related to the locally induced gap = —U {d^^d-^)). 
They read 


U U ^ 

Ed-e+---J2 


E, 
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i5 = r cos 


$ 


--y. 


^ — r cos Y 


1 - 


A 




D^^{iuJri) 


(9) 


Since we are primarily interested in the zero-temperature 
QPT where the energies of the ABS approach zero ujq —>■ 
0, we can approximate the denominators in the integrals 
by their low-frequency asymptotics 0) Ki E‘^ + 

(5^ -I- (1 -I- r/A)^a;^, which implies wq ~ 4- (5^/(1 -I- 

r/A). Near the quantum critical point we then obtain 


Ed 

5 


1-f 

1-f 


U 


20.0 (1 + i ) 2 

U 

2^0 (1 + ^ 


U 


= T cos ■ 


4) 




( 10 ) 


with the band contribution expressed via the function 

7-/™ (fi'i _ dt _ cosh^ t _ 

1 . Jo 27r cosh^(i/2)(a:-t-cosh t)^+x^ cos^(<&/2) sinh^(i/2) ' 

Re-parametrizing 

Ed = (1 + r /A)uJo cos 4’ and d = (1 + F / A)ujq sin 4 we 
arrive at 


U 


2(l + £)^ 

U 

2(1 + D' 


■ 0.0 


- OJQ 


u 


/ 2 
COS 4 = 


1-f 


r ’ 


r cos „ 

sin'0 = 


<s> r 


1-f 


1-f 


A^ 




( 11 ) 


At the QPT characterized by 0.0 = 0 the solubility 
condition (cos^ 4 + sin^ ^/> = 1) gives us the equation for 
the HF phase boundary 


u 

2 

r u] 

.2(l + i). 


2 J 


+r^ cos^ 


$ 




$ 


1 2 


( 12 ) 

that generalizes the corresponding well-known expression 
in the atomic limit A —>■ 00 [ 23 , [13 • This HF phase 
boundary plotted in Fig. [2] is not particularly precise, 
however, it yields qualitatively reasonable results. More¬ 
over, we have noticed that when the band contribution I 
in equation (1121) is omitted one gets a surprisingly good 
and extremely simple approximation for the boundary, 
that we call here the generalized atomic limit (GAL), ly¬ 
ing for half-filling (e = —U12) typically very close to the 
numerically exact results by NRG, see Fig. [2^-b. Obvi¬ 
ously, the HF approximation heavily overestimates the 
contribution from the band continuum. 

Eqs. (fTTI) may be used also around the QPT, when ujq 
is small (and unknown). We can see that wq is positive 
on one side of the boundary while it is negative on the 
other side. Since wq > 0 by construction, we must con¬ 
clude that the solution with negative ojq, that we iden¬ 
tify with the TT-phase region, is unphysical. We cannot 
go beyond the phase boundary from the 0-phase to the 
TT-phase within this perturbative approach based on the 
assumption of a nondegenerate ground state. 
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Figure 3. Comparing various methods of calculation of one- 
particle quantities. Panels a and b show supercurrent at half¬ 
filling as a function of the phase difference <f> for U = 4A (a) 
and U = 8A (b) calculated by numerically exact NRG, and 
analytically approximative fRG, spin-symmetric HF and, fi¬ 
nally, the second-order PT/dynamical corrections (DC) show¬ 
ing a nearly perfect agreement with NRG (unlike the other 
two methods). Inset in panel a depicts the ABS energies ojo as 
functions of <E> for the two values of the Coulomb interaction 
U. The green dashed line in panel b represents the HF tun¬ 
neling current component. In panel c the occupation number 
n = (dlrdtr) and locally induced SC gap Ad = —U {d_id4 (in¬ 
set) are plotted as functions of the level energy for two values 
of the phase difference 4> = 0 (with no phase transition) and 
>1> = TT (exhibiting phase transition). fRG data in panels a, b 
were graphically taken from Fig. 4b of Ref. f^ . 
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Dynamical corrections 


The qualitative predictions of the HF approximation 
can be significantly improved by including dynamical cor¬ 
rections into the self-energy, which come from the second 
order of the perturbation expansion represented by the 
second and third diagrams on the r.h.s. of equations in 
Fig-dl The two diagrams originate in two different types 
of intermediate propagation consisting of either normal 
or anomalous propagator. The mathematical equivalents 
for the second-order contributions read 

= ^ GiiUJn + (13) 

m£Z 

and 

Jj2 

= J- ^ -f Wm)xii^m) (14) 

m£Z 


where 

XiiVm) 


1 


[G(iuJn)G{iU}n + iVm) + G {iU!n)G(iuJn + Wm 

n^Z 


(15) 


is the two-particle bubble consisting of the normal and 
anomalous parts and Vm = (3 is the m-th bosonic 

Matsubara frequency. 

These first two orders of the perturbation expansion 
are well controllable on the one-particle level. The higher 
contributions to the self-energy become more complex 
and their classification more complicated. For a general 
discussion of this problem see Ref. [sij . 

The second order self-energy correction together with 
the first-order (in U) HF counterparts are inserted into 
the equation for the Green function, equation We 
obtain a self-consistent nonlinear functional equation for 
the whole Green function as a function of frequency. This 
equation is solved numerically at zero temperature. We 
noticed, however, that nearly identical results are ob¬ 
tained by computationally less elaborate method which 
evaluates the dynamical self-energies by using just a fully 
converged HF solution as the input GF. The convolutions 
in the second-order self-energies are thus evaluated just 
once at the beginning of the procedure and consequently 
used as fixed inputs into the self-consistent procedure it¬ 
erating the Green function through the HF self-energy. It 
should be stressed that while the second-order contribu¬ 
tion may be simplified in this way, the full self-consistency 
loop between the GF and the HF self-energy is manda¬ 
tory — any compromises there lead to even qualitatively 
wrong results. 


DISCUSSION 

We have carried out the above mentioned procedure 
both in the Matsubara formalism as well as in the spec¬ 
tral representation (performing the analytic continuation 


described in Methods) with identical results. We have 
found that the 0-phase smoothly develops from the non¬ 
interacting limit U = 0 and terminates at the 0 — tt 
phase boundary beyond which there exists no regular self- 
consistent solution for the GF. In the spectral represen¬ 
tation this is associated with the energy of ABS ojq reach¬ 
ing zero. The results for the phase boundaries, shown in 
Fig. [21 and one-particle quantities in the 0-phase in Fig. [3] 
exhibit unprecedented precision of the dynamical correc¬ 
tions approximation which gives numerical results nearly 
identical to the numerically exact NRG data produced 
by the “NRG Ljubljana” open source code for all 

studied parameter sets as well as all physical quantities. 
Surprisingly, it outperforms in the regime of not-so-weak 
interaction even the fRG method designed for the strong 
correlations (see the H-axis scale in Fig.|2^). This is likely 
due to the static-vertex implementation of the fRG in 
Ref. (13 • The limitations of the static-vertex approxima¬ 
tion have been discussed before (see Ref. [s^, Sec. 9.4.6), 
nevertheless it is currently the only one technically viable 
for fRG. On the other hand our dynamical corrections 
)]properly include the frequency dependence of the corre¬ 
lation effects (even if just perturbatively) which probably 
explains their superiority over the fRG in the description 
of 0-phase quantities as well as the phase boundary. In 
this context, we would also like to point out an interest¬ 
ing observation we have made. In Fig. |3]3 we plot (by the 
green dashed line) the tunnelling part of the supercurrent 
(the second term in the lower equation (|7])) for the HF 
solution and see that it coincides in the overlapping range 
of parameters with the full supercurrent solution of the 
fRG in the 7r-phase. Although plotted for clarity just in 
Fig.lSb this observation holds for all J —$ characteristics 
taken graphically from Ref. Since our HF solution 
breaks down at the phase boundary we cannot extrapo¬ 
late beyond it, nevertheless there is obviously some subtle 
correspondence between the spin-symmetric HF solution 
and the 7r-phase solution of the fRG. 

To conclude, we have presented a systematic perturba¬ 
tive expansion for the 0 — tt transition in the supercon¬ 
ducting Anderson model and found out that its second 
order yields at zero temperature excellent results for the 
phase boundary and quantities in the 0-phase such as lo¬ 
cally induced superconducting gap or supercurrent sur¬ 
passing any (semi) analytical methods employed to this 
model so far. Although demonstrated here explicitly 
just for the symmetric case T^ = F/j for simplicity, the 
method produces equally good results also in the general 
case. Moreover, we have also verified numerically that 
the formalism is gauge-invariant, i.e., physical quantities 
depend on the phase difference — <!>/{ only and con¬ 
serves current, i.e., supercurrents calculated at left/right 
junctions are identical. Furthermore, the full second- 
order PT is thermodynamically consistent (unlike, e.g., 
fRG H). 

The method cannot be, however, continued to the tt- 
phase without modifications taking into account the de¬ 
generacy of the doublet ground state. Moreover, we 
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have observed that the Matsubara formalism at finite 
temperatures does not detect any sharp phase bound¬ 
ary found at T = 0. To our best knowledge there is 
presently no (semi) analytical method that would concep¬ 
tually correctly and quantitatively reasonably describe 
the TT-phase. The spin-polarized HF suffers from the dis¬ 
continuity problems mentioned in the Introduction while 
the fRG solution returns e- and U-independent quanti¬ 
ties in the 7r-phase apparently closely related with 
the simplest spin-symmetric Hartree-Fock approximation 
as discussed above, which is clearly not sufficient. The 
construction of an analytical theory of the 7r-phase thus 
remains an open challenge for future study. 


METHODS 

The necessary information for the study of the crossing 
of ABS as well as for obtaining the particular components 
of total current can not be obtained directly from the ex¬ 
pressions in Matsubara frequencies. To access it we ana¬ 
lytically continued the expressions to the real-frequency 
domain. 

The inverse Green function ([4]) can be represented as 




- 4 - 3 - 2-101234 
co/A 


- 1 /z[l-fs(z)]-e-S(z) A$s(z)-5(z) \ 

V A$s(z) —5(—z) z[l -I- s(z)] -I- e -I- 1](—z)/ 

(16) 

where 


s(z) 


iV 


sgn(Im z) 


(17) 


is a dynamical renormalization of the impurity energy 
level due to the hybridization to the superconducting 
leads. We introduced a renormalized complex energy 
C = ^ + related to z = a; -I- zy via = z^ — A^. The 
following convention for complex square root is used: 


^r] = ojy, sgn^ = sgnw, sgnr] = sgny, (18) 

so that C = z for A = 0. The renormalized energy C 
along the real axis is then real outside the energy gap 
and imaginary within it. Accordingly to this definition 
the function s(z) is imaginary outside the energy gap and 
real within it, 


s{uj ± zO) 

zF sgn(a;) 

for 

<1 

A 

\/w2 - A2 



s{ljJ ± zO) 

i 

for 

<1 

V 

VA^ — 




This definition allows for a straightforward analytic con¬ 
tinuation of the Matsubara Green function to real fre¬ 
quencies. An illustrating example of the normal and 
anomalous spectral functions is plotted in Fig. ID The 
Green function has a gap around the Fermi energy from 
—A to A and two poles at iwo, IwqI < A. The positions 
of these poles are given by zeroes of the determinant. 


Figure 4. Normal (— Im G/tt, upper panel (a)) and anomalous 
{—ImO/Tv, lower panel (b)) spectral density for U = 4A, 
r = 2A, <1> = 7r/2 and e = —V12 (half-filling) calculated 
using the dynamical corrections from the second-order of the 
perturbation expansion. The heights of the arrows marking 
the Andreev bound states represent their residues. 


det[G'“^(tLio)] = 0. Since the function s(uj) has a square- 
root singularity at gap edges, the gap is fixed and does 
not depend on interaction strength. 

Galculating the self-energy from diagrammatic expan¬ 
sion calls for the analytic continuation of sums over Mat¬ 
subara frequencies. The sum of a one-particle function 
F over fermionic Matsubara frequencies can be rewritten 
in the spectral representation as 

dujfiuj) ImF(w -I- zO) 

( 20 ) 

where uji are the isolated poles within the gap and 
/(w) is the Fermi-Dirac function. This formula can be 
used directly to calculate the static Hartree-Fock self¬ 
energies dSJ and the Josephson current d?]). 

Similar approach can be utilized to calculate the two- 
particle bubbles and the second-order dynamic correc¬ 
tions, Eqs. (fni)-([TS|). For the sake of simplicity we resort 
to zero temperature. Ghoosing a correct contour in the 
upper complex half-plane we arrive at an expression for 


1 —A . , 1 

- F{lUJn) ^ - 


-A 


/A 
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the normal part of the bubble, 


1 

Xn{ijJ~^) = - J dxlm.G{x'^)\G{x + w'^) + G{x — uj'^)\ 

+ Res(G', —wo)[C (—wq “f C!{—ujq — w)] 

( 21 ) 


and analogously for the anomalous part Xa- We have 
abbreviated = u! + iO. The resulting bubble has an 
extended gap from — A — wq to A + wq- The contribu¬ 
tions from the isolated states at ±2wo from the normal 
and anomalous parts exactly cancel out each other, so 
there are no gap states in the full bubble y = Xn + Xa- 
Taking this into consideration we arrive at a formula for 


the normal self-energy, 

= — f dxlmG{x^)x{x — oj^) 

J — OO 

tj 2 p-A-u}Q ^22) 

H-/ dx Im xix^)G{x + uj~^) 

- Res(G', -a;o)x(-wo - w) 

and similarly for . Integrals of this kind can be evalu¬ 
ated numerically using fast Fourier transform algorithms 
which makes the calculation simple and efficient. 
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